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We investigate ground state and finite temperature properties of the half-filled Hubbard model on 
a honeycomb lattice using quantum monte carlo and series expansion techniques. Unlike the square 
lattice, for which magnetic order exists at T — for any non-zero U, the honeycomb lattice is known 
to have a semi-metal phase at small U and an antiferromagnetic one at large U. We investigate the 
phase transition at T = by studying the magnetic structure factor and compressibility using quan- 
tum monte carlo simulations and by calculating the sublattice magnetization, uniform susceptibility, 
spin-wave and single hole dispersion using series expansions around the ordered phase. Our results 
are consistent with a single continuous transition at U c /t in the range 4 — 5. Finite temperature 
signatures of this phase transition are seen in the behavior of the specific heat, C(T), which changes 
from a two-peaked structure for U > U c to a one-peaked structure for U < U c . Furthermore, the U 
dependence of the low temperature coefficient of C(T) exhibits an anomaly at U ~ U c . 

PACS numbers: 71.10.Fd,75.10.Lp,75.40.Mg 
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I. INTRODUCTION 

The two-dimensional Hubbard Hamiltonian has been 
extensively studied as a model of metal-insulator and 
magnetic phase transitions [J and also within the con- 
text of systems such as the CuC>2 sheets of high tem- 
perature superconductors. [2 In the square-lattice case, 
at half-filling, nesting of the Fermi surface leads to a di- 
vergent antiferromagnetic susceptibility as the tempera- 
ture is lowered, even for U = 0, and thus the ground 
state is an antiferromagnetic insulator at any non-zero 
U. It is of interest to study cases where, instead, the 
transition to the antiferromagnetic phase occurs at finite 
U. In such a situation, for example, it may prove possi- 
ble to see if the Mott metal- insulator and paramagnetic- 
antifcrromagnetic phase transitions occur separately. A 
finite U c also makes it more straightforward to study the 
thermodynamics at temperatures above the T — quan- 
tum phase transition, a question pertinent to experimen- 
tal studies of such phase transitions. 

The two-dimensional honeycomb lattice is one geom- 
etry in which we can explore these issues. In this pa- 
per, we investigate the properties of the half-filled hon- 
eycomb lattice Hubbard model using determinant quan- 
tum monte carlo and series expansions methods. After a 
brief review of previous work, we describe the model and 
calculational approaches, and show data for a number 
of different ground state properties that carry signatures 
of the phase transition. Our overall results are consis- 
tent with a single continuous transition as a function of 
U/t. We then turn to the finite temperature behavior of 
the specific heat to see how such a critical point may be 
reflected in this key experimental property. 

While the honeycomb lattice has U c non-zero, it is un- 



its non-interacting density of states has a special feature. 
As shown in Fig. 1, N(lu) vanishes linearly as u> — > 0, 
so the system is a semi-metal (or alternatively, a zero- 
gap semiconductor) at half-filling. As a consequence, at 
weak coupling, the low temperature behavior of the spe- 
cific heat is quadratic in temperature, C — ST 2 , instead 
of the usual linear Fermi liquid dependence. At strong 
coupling, when long range antiferromagnetic order sets 
in, the specific heat will also be quadratic in T owing to 
the spin-wave excitations. How the specific heat evolves 
between these two regimes is an open question. 

A considerable body of work exists concerning the 
ground state phase diagram. Martelo et al found that 
within mean field theory the Mott transition occurs be- 
low an upper bound for the critical interaction strength 
U c /t fa 5.3. |H Meanwhile, their variational monte carlo 
calculation suggested a lower bound for the antiferomag- 
netic transition U c /t ~ 3.7. They interpreted these re- 
sults as a single transition from paramagnetic metal to 
antiferromagnetic insulator at U c /t = 4.5 ± 0.5. 

Baskaran et al |j| and Sorella et al J5| studied the 
model using the Random Phase Approximation which 
gives U c /t = 2.23 for the onset of antiferromagnetic 
order. Associated auxiliary field quantum monte carlo 
(QMC) simulations 3 using the ground state projec- 
tion approach suggested U c /t — 4.5 ± 0.5. Later QMC 
work by Furukawa[6j on larger lattices and doing system 
size extrapolations resulted in a somewhat lower value, 
U c /t = 3.6 ± 0.1. Peres et al have recently studied the 
phase diagram and mean field magnetization of coupled 
honeycomb layers as a function of filling, U/t, and inter- 
layer hopping t' /t using the random phase approximation 
and spin wave analysis. 7] 

As with the square lattice Hubbard model, Nagaoka 




FIG. 1: The noninteracting density of states of the Hubbard 
model on a honeycomb lattice. This geometry is bipartite, so 
N(u>) is symmetric about u> = 0. The vertical lines correspond 
to fillings of 0.1, 0.2, 0.3, . . . The density of states vanishes 
linearly at u> — and has logarithmic Van Hove singularities 
at fillings p = 3/8 and 5/8. The bandwidth W = 6t 



lattice at strong couplings (U/t > 12), as has been inves- 
tigated by Hirsch using exact diagonalization on small 
clusters |3 and by Hanisch et al using Gutzwiller wave 
functions. J3] 

We conclude this introduction by noting that the Hub- 
bard model on a honeycomb lattice has also been sug- 
gested to be of interest for a variety of systems including 
graphite sheets, Q and carbon nanotubesjlCJ]. as well as 
MgB 2 (3 and Pb and Sn on Ge(lll) surfaces. |ll| The 
honeycomb lattice is also a 2/3 subset of the triangu- 
lar lattice, and so the nature of spin correlations on the 
honeycomb lattice has been considered as possibly rele- 
vant to the properties of triangular lattice systems like 
Na^CoG^ at appropriate dopings |12 |. 



II. THE HUBBARD HAMILTONIAN, 

DETERMINANT QUANTUM MONTE CARLO, 

AND SERIES EXPANSION METHODS 

We study the Hubbard Hamiltonian, 



H = ^E^t^ + i^) 



<ij><7 



uJ2( n n - ^)(«u - 5) -m^(«it + n u)- 



Here c[ a (c-) are creation(destruction) operators for a 
fermion of spin a on lattice site i. The kinetic energy 
term includes a sum over near neighbors (i, j) on a two- 
dimensional honeycomb lattice. We denote by N the 
number of lattice sites, and L the linear dimension with 
N = |L 2 . The interaction term is written in particle- 



filling: the density p = (n^ H-nij.) = 1 for all t, U and tem- 
peratures T. We choose the hopping parameter t = 1 to 
set the energy scale. Note that the noninteracting model 
has two 'Dirac' points K± on the Fermi surface where 
the dispersion relation is relativistic,jj| i.e. the energy 
grows linearly with | k — K± | . 

We use the determinant quantum monte carlo method 
to measure the properties of the Hamiltonian. [Tj|] In this 
approach the partition function is written as a path in- 
tegral, the interaction term is decoupled through the in- 
troduction of a discrete auxiliary 'Hubbard-Stratonovich' 
field, [lj] and the fermion degrees of freedom are traced 
out analytically. The remaining summation over the 
Hubbard-Stratonovich field is done stochastically. Since 
the lattice is bipartite, no sign problem occurs at half- 
filling. Data were typically generated by doing several 
tens of thousands of measurements at each data point 
(temperature, coupling constant, lattice size). 'Global 
moves' which flip the Hubbard-Stratonovich variables for 
all imaginary times at a given spatial site were included 
so that at stronger couplings, transitions between differ- 
ent densities are facilitated. 15] 

We have also carried out an Ising type expansion for 
this system at T = using a linked-cluster method. [la] 
Similar expansions were previously done for the Hubbard 
model on the square lattice. [ItJ To perform the series ex- 
pansion, one needs to introduce an Ising interaction into 
the Hubbard Hamiltonian, and divide the Hamiltonian 
into an unperturbed Hamiltonian (Ho) and a perturba- 
tion (H\) as follows, 



H = H + XHi 

Ho 



J E( CT i °l + J ) + ^EK - h^i b 



ff l = X> JKoj* + !) - *(4r c J* + h ' C -)] 
(«> 

where of = m-\ — nij., and A is the expansion parameter. 
Note that we are primarily interested in the behavior 
of the system at A = 1, at which point the Ising term 
cancels between Ho and H\ . The strength of the Ising 
interaction can be varied to improve convergence, and 
it proves useful to keep it of order t 2 /t/.[T3 The limits 
A = and A = 1 correspond to the Ising model and the 
original model, respectively. The unperturbed ground 
state is the usual Neel state. The Ising series have been 
calculated to order A 15 for the ground state energy, the 
staggered magnetization M, and the square of the local 
moment I>m[ljj], an d to order A 13 for the uniform mag- 
netic susceptibility \±- The resulting power series in A 
for t/U = 0.15 and J/U = 0.0225 are presented in Table 
I. 

In addition to the ground state properties, we also com- 
pute the spin- wave dispersion A (to order A 13 ) and the 
dispersion An of 1-hole doped to half-filled system (to 
order A 11 ). The calculation of the dispersion involves a 
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FIG. 2: Density p as a function of chemical potential p at 
weak coupling (U — 2t). p is not pinned at one, but immedi- 
ately begins to shift when p 7^ 0: there is no Mott gap. 



FIG. 3: Density p as a function of chemical potential p at 
strong coupling (U — 7t). There is clear evidence for a Mott 
gap. 



dispersions can be written in the following form 
A(k x ,k y ) = ^2^2-aij^X p {cos(ik x /2)cos(V3jk y /2) 

+ cos[(i + 3j)fc x /4] cos[V3(i - j)k y )/4] 
+ cos[(i - 3j)fc x /4] cos[V3(i + j)k y /A] } 

In Table II, we list the series coefficients a,.^ for t/U = 
0.15 and J/U = 0.0225. The series for other couplings 
are available from the authors upon request. 

III. QUANTUM PHASE TRANSITION 

We begin by examining the evidence for a phase tran- 
sition in the model. First, we present results from the 
quantum monte carlo simulations, which can, in princi- 
ple, address arbitrary t/U ratios. 

A. Compressibility 

The Mott metal-insulator transition is signalled by a 
vanishing compressibility k = dp/ dp. We show p as a 
function of p for (3 = 8 and U = 2* (Fig. 2) and U = It 
(Fig. 3). There is a clear qualitative difference in behav- 
ior. For weak coupling, p immediately shifts from half- 
filling as p increases from zero, while at strong coupling, 
p remains pinned at p = 1 out to finite p. 

In Fig. 4 we show the filling as a function of U at a 
small non-zero value of the chemical potential po = 0.20. 
We see that po — > 1 at U « bt, signalling the onset of the 
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FIG. 4: The difference in the value of the density p from half- 
filling at small nonzero po = 0.2i is a measure of the presence 
the Mott gap. Here we see p(po) — » 1 and hence a gap opens 
at U/t « 5. 



B. Spin Correlations and Antiferromagnetic 
Susceptibility 

To study the magnetic behavior, we measure the real 
space spin correlations, 

c**(r) = (S z (r)S z (0)) S z (r) = n rT - n H 

c+_(r) = (5_(r)5+(0)) S+(r) = c^c ri , 

and their Fourier transforms, 
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FIG. 5: The antiferromagnetic structure factor is shown as a 
function of inverse temperature /3 and different lattices sizes 
L at U/t = 7. At low f3 (high T) the correlation length is 
short, and S(tv,tt) is independent of L. At large /3, S(tv,tt) 
grows with L. 



At T = and in the antiferromagnetically ordered 
phase at large U/t, the real space correlation will go 
asymptotically to a non-zero value m 2 /3 at large separa- 
tions r. In our finite temperature simulations, we access 
the T = limit by cooling the system to the point where 
the correlation length exceeds the lattice size. In this 
case, the structure factor will grow linearly with lattice 
size N. More precisely, the structure factor will obey, 



TV 



S(q) =m 2 /3 + a/L, 



where L is the linear lattice size. |lj| In the paramagnetic 
phase at small U/t, the structure factor will be indepen- 
dent of N , and hence S(q)/N will vanish as N — * oo. 

In Fig. 5 we show S(tt, it) as a function of inverse tem- 
perature /3 for different lattice sizes at U/t = 7. Fig. 6 
shows the associated scaling plot. Also shown in Fig. 6 
is the value of c(r) for the largest separations on our fi- 
nite lattices. This quantity should scale with the same 
intercept m 2 /3 but a different finite size correction. We 
see that the system is in an antiferromagnetically ordered 
phase for this coupling. 

Figs. 7-8 show analogous plots at U/t = 6. The order 
parameter is still non-zero, but is quite small. Similar 
plots for U/t = 5 are consistent with the vanishing of 
long range order. While we cannot pin down the location 
of the quantum phase transition exactly, a comparison of 
this analysis with the compressibility of Fig. 3 suggests 
that the vanishing of the compressibility gap and the an- 
tiferromagnetic order occur very close to each other and 
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FIG. 6: The scaled structure factor (filled triangles) and 
spin correlation function (filled squares) at large distance are 
shown for large j3 as a function of the inverse linear dimen- 
sion for U/t — 7. The lines are least squares fits to the data. 
These quantities scale to a nonzero value of the order param- 
eter (square of the staggered magnetization) in the thermo- 
dynamic limit 1/L — > 0. 
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FIG. 7: Same as Fig. 5 except U/t = 6. 

C. Results from Series Expansions 

We now present results from the Ising type series ex- 
pansions. These expansions are only valid in the magnet- 
ically ordered phase, and thus can only access the prop- 
erties of the system for U > U c . 

In Fig. |5| and Fig. ^] we show the sublattice magneti- 
zation and uniform susceptibility. The QMC results for 
the sublattice magnetization are also shown. The two 
agree with each other for small t/U. The uncertainties 
increase as the transition is approached. QMC results 
suggest a more abrupt drop to zero around U/t ss 5, 
whereas the series results suggest a gradual decrease with 



d 



oo d 



_ E d 



o 



A - 

▲ \/^ 
■ ' ■ ^^ 
A /^ 







0.1 
1/L 



0.2 



FIG. 8: Same as Fig. 6 except U/t — 6. 




FIG. 9: The staggered magnetization versus t/U obtained 
from series expansions and quantum monte carlo simulations. 
The lines joining the points are a guide to the eye. See text 
for more discussion. 



variable t/U but rather in an auxiliary variable A, it is 
difficult to locate the true critical point U c /t and obtain 
the critical properties. However, since we expect the crit- 
ical exponent /3 to be less than one, the true curve should 
come to zero with an infinite slope. Thus, from the series 
results alone, one would estimate U c /t sa 4, and this is in 
agreement with the estimate from the susceptibility x± 
shown in Fig. I1UI 

Next, in Fig. 1111 we show the spin-wave dispersion 
along high-symmetry cuts through the Brillouin zone for 
t/U = 0.1, together with the dispersion obtained from 
first and second order spin-wave results for the Heisen- 
berg model on a honeycomb lattice[ljj, which should ap- 
proach the dispersion for the Hubbard model in the large 
U limit. We can see that the dispersion has its minimum 




FIG. 10: The uniform susceptibility U\± versus t/U obtained 
from series expansion. 



berg model on a honeycomb lattice has a maximum at 
W point, while for the Hubbard model, this is only true 
for very small t/U. Already for t/U — 0.1, the energy at 
W points is reduced, and the maximum moves to the K 
point. 

Also, in figure 1121 we show the 1-hole dispersion for 
selected values of t/U, where we can see that the min- 
imum and maximum gaps are at the W and T points, 
respectively. This dispersion is quite different from the 
case of the square lattice, since there is no nesting of the 
Fermi surface here. For the square lattice, the single hole 
dispersion relation is anomalously flat near the degener- 
ate points (0, ±7r), (±7t, 0) of the Brillouin zone, with the 
minimum of the dispersion at (±7r/2,±7r/2).[2£J Fig. IT31 
shows the minimum gap, i.e. the gap at the W point, 
and the bandwidth, Ar — A^y, vs t/U. The gap closes 
at t/U ~ 0.26, indicating a transition to the semi-metal 
phase. 

To summarize, study of both magnetic and charge 
properties using series expansions show a direct transi- 
tion from the antifcrromagnetic to the semi-metal phase 
around U c « At. 

Combining the quantum monte carlo and series expan- 
sion results, we estimate the phase transition to be in the 
range U c /t = 4 — 5. There is greater internal consistency 
in the location of the critical point if we restrict ourselves 
to one method. But, in fact, there are larger uncertainties 
in both methods especially as the quantum phase transi- 
tion is reached. However, both methods strongly indicate 
that the Mott transition and the antiferromagnetic order 
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FIG. 11: Plot of the spin-wave excitation spectrum A(k x , k y ) 
(in units of effective J e g = 4i 2 /(7) along the path YKWY in 
the Brillouin zone (see the inset, where the momentum k for 
T, K, W points are (0,0), (2tt/3,0), and (2tt/3, 2tt/3>/3), 
respectively) for the system with coupling ratios t/U = 0.1 
in the Neel ordered phase. Also shown are the first (das hed 
line) and second (solid line) order spin- wave results [l9|| for 
Heisenberg model on honeycomb lattice. 
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FIG. 12: Plot of the 1-hole excitation spectrum A(k x ,k y )/U 
in the Neel ordered phase along the path TKWY in the Bril- 
louin zone (see the inset) for the system with coupling ratios 
i/[/ = 0.05, 0.01, 0.15, 0.25. 



FIG. 13: The minimum single-hole gap Aw at W point and 
its bandwidth Ar — Aw vs t/U. 



IV. SIGNATURES OF THE QUANTUM PHASE 
TRANSITION IN THE SPECIFIC HEAT 



An important objective of our study was to examine 
the signature of the quantum phase transition in the fi- 
nite temperature behavior of the specific heat. We now 
turn to those studies, which are based on the quantum 
monte carlo method. 

At strong couplings, one expects two features in the 
specific heat of the Hubbard Hamiltonian. The first, at a 
temperature T « U/5, signals the formation of magnetic 
moments. [2ll |22| while the second, at a lower tempera- 
ture T s=s J = At 2 /U, is associated with the entropy of 
moment ordering. This picture has been verified in the 
one-dimensional case using Bethe Ansatz techniques [23| 
and (u sing quantum monte carlo) in the two dimensional 
square [24 l25l | and three dimensional cubic lattices. [261] 
Interestingly, in the square lattice, the two peak struc- 
ture persists to weak coup ling where the energy scales 
U and J have merged. J2l| In one dimension, there is a 



single peak at weak coupling. [27ll28| 

The specific heat C(T) for the two-dimensional honey- 
comb lattice is shown in Fig. 16 for different couplings U 
and lattice size L — 12. For strong coupling, U/t = 6, 7, 8 
there is a clear two peak structure. This is replaced by 
a single peak for weaker couplings, U/t = 2, 4, 5. Again, 
this result is in contrast with the behavior of C(T) on the 
square lattice, where a two-peak structure is evident for 
all C//T.[2l| It is plausible to conjecture that the differ- 
ence is the absence of long range antiferromagnetic order. 
This suggestion is supported by the fact that coalescence 
of the specific heat peaks is seen in "Dynamical Mean 
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FIG. 14: The specific heat C(T) is shown as a function of 
temperature for different coupling strengths. In the antiferro- 
magnetic phase for U > U c , the specific heat has a two peak 
structure. In the metallic phase for U < U c there is a single 
peak. The 'universal crossing' at T = 1.6£ is discussed in the 
text. 



are restricted to the paramagnetic phase and antiferro- 
magnetic fluctuations are neglected. 

This is, however, a rather subtle question, since the 
Mermin- Wagner theorem precludes long range order at 
finite temperature. What is meant, more precisely, is 
that on a two-dimensional square lattice, the low T struc- 
ture in C(T) appears when the antiferromagnetic correla- 
tion length £(T), begins to grow exponentially as T — ■> 0. 

The evolution from a two to a one peak structure 
in C(T) is one interesting reflection of the underlying 
quantum phase transition on the finite temperature ther- 
modynamics. Another way of examining this question 
concerns the low temperature behavior of C(T). As 
pointed out in the introduction, we expect a quadratic 
temperature dependence at both strong coupling (spin- 
waves) and weak coupling (linearly vanishing density of 
states at the Fermi level). How does the coefficient 8 
in C(T) = 8T 2 evolve as one crosses between the two 
phases? 

Before we present the results for 8, we note that ex- 
tracting 5 is clearly a subtle numerical issue. On the 
one hand, 8 characterizes the low T behavior, but on the 
other hand, because of finite size effects, which become 
larger as the temperature is lowered, one cannot use data 
at too low values of T. Thus, our calculation of 8 should 
be viewed with some caution. What we have done in gen- 
erating Figs. 15-16 is to fit the data for C(T) to the T 2 
form over only a finite temperature window: below the 
peak in C(T) but also above the temperatures at which 
finite size effects begin introducing a noticeable gap in 
the spectrum. In Fig. 15 we show 8 as a function of U. 



value U c /t sa 5 previously inferred from the compress- 
ibility and spin correlation data. Fig. 16 emphasizes this 
feature by plotting the derivative of 8 with respect to 
t/U as a function of t/U. As we have noted, the specific 
heat of the noninteracting system obeys C = S(U = 0)T 2 
with 8(U = 0) = 4.1, because of the linearly vanishing 
density of states. Perturbation theory suggests that for 
small finite U, 8 should increase quadratically from this 
value. Nevertheless, in the vicinity below the quantum 
phase transition, the value for 8 extracted from the quan- 
tum monte carlo data looks rather linear in U, as seen in 
Fig. 15. If 8 = mU/t then d{8)/d{t/U) = -m/(t/U) 2 . 
With this in mind, a line showing the functional form 
—m/(t/U) 2 with m = 2 is given and fits the weak cou- 
pling data very well. The breakaway from this form at 
strong coupling further emphasizes the change in behav- 
ior in the vicinity of the quantum phase transition. 

In studies of the two peak structure of the specific heat 
on the square lattice, an interesting interchange of the 
role of kinetic and potential energies was noted. [21j At 
large U, the temperature derivative of the potential en- 
ergy was the primary contribution to the high T, 'mo- 
ment formation', peak, while the temperature derivative 
of the kinetic energy drove the low T, 'moment ordering', 
peak. However, at weak U the situation was reversed, 
with the high T peak originating in the kinetic energy. 
With that separation in mind, we plot in Fig. 17, for 
the honeycomb lattice, the contributions of the potential 
and kinetic energies to 8. It is the contribution of the 
potential energy to 8 which appears to have the sharper 
evolution in the vicinity of the quantum phase transition. 

Returning to the specific heat versus temperature, 
shown in Fig. 14, we note the existence of a very well 
defined crossing point at T ss l.6t. This crossing has 
been observed previously in DMF T.I29. J3(l l31| and in 



the two dimensional square lattice. [2 lll24| Indeed, in the 
former case, two crossings were observed, with the high 
temperature one being nearly universal, while the low 
temperature intersections were considerably more spread 
out, much as we observe in Fig. 14. It is also interesting 
that the numerical value of the crossing is almost identi- 
cal for the honeycomb and square lattices, despite their 
different bandwidths. 

Finally, we turn to the behavior of the entropy 5*. In 
Fig. 18 we show S as a function of U for different tem- 
peratures T. At large U, the clustering of the curves 
for different temperatures near ln(2) is indicative of the 
existence of disordered magnetic moments in a range of 
intermediate T . The low temperature magnetic ordering 
tendency is evident in the gap between the T = 0.2 and 
T = 0.3 curves. As U is decreased, the screening away of 
the moments is indicated by the T = 0.3 isotherm drop- 
ping from ln(2) to 0. It is interesting that this behavior 
is so gradual. Finally at small U one observes the more 
or less equally spaced isotherms of free electron gas. This 
figure complements the data of C(T) shown in Fig. 14, 
since the entropy hang up at large U near ln(2) is just 
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FIG. 15: 8, the coefficient of the T 2 term in the specific heat 
is shown as a function of U/t. The solid square is the U — 
value. There appears to be a change in slope as U crosses U c . 



FIG. 17: The separate contributions of the potential (6u) 
and kinetic (<5jf ) energies to the quadratic coefficient of the 
specific heat are shown. Su shows the more abrupt behavior 
in the vicinity of U c . The small differences between the values 
of 5 obtained from the total energy, and the values Sk + 8u 
from the kinetic and potential energies separately provide a 
measure of the uncertainties in our fitting procedure. 
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FIG. 16: The derivative of S, the coefficient of the T term in 
the specific heat, with respect to t/U is shown. This derivative 
has a sharp change near the critical coupling U c - The solid 
line is -2/(t/U) 2 (see text). 



FIG. 18: The entropy is shown as a function of U for different 
temperatures. At large U the gaps between the T — 10 and 
T — 2 curves and between the T = 0.3 and T = 0.2 curves 
reflect the entropy loss associated with magnetic moment for- 
mation and ordering respectively. 



Figure 19 exhibits the entropy as a function of tem- 
perature. At weak coupling, there is a smooth evolution 
from ln(4) at high T to zero at low T. For strong cou- 
pling, a plateau near ln(2) interrupts this evolution, again 
exhibiting a range of temperatures with well formed, but 



V. CONCLUSIONS 

In this paper we have studied the Hubbard Hamilto- 
nian on a half-filled honeycomb lattice using quantum 
monte carlo and series expansion methods. Both meth- 
ods strongly suggest that the model has a single con- 
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FIG. 19: The entropy is shown as a function of T for weak 
and strong coupling. The dashed and solid lines are the results 
of quantum monte carlo simulations on 12x12 lattices. The 
symbols are generated by an exact calculation on a two site 
model for comparison. 



finite temperature signatures of the phase transition in 
the specific heat, as a guide to experimental studies. We 
observe that around U c the specific heat changes from 
a one peak (below U c ) to a two peak (above U c ) struc- 
ture. We suggest that this is associated with the fact 
that for U > U c the antiferromagnetic correlation length 
grows rapidly as the temperature is reduced. For weak 
coupling only very short-range antiferromagnetic corre- 
lations exist, and the specific heat has no signature of 
magnetic order. 

We also studied the evolution with on-site interaction 
strength U of the coefficient 6(U) of the quadratic tem- 
perature dependence of the specific heat at low temper- 
atures. Since the excitations which produce the T 2 term 
above and below the quantum phase transition are un- 
related, one might have expected S(U) to exhibit a dis- 
continuity at U c . Instead, we found a sharp change in 
the slope, d5(U)/dU at U c . Given the uncertainties in 
obtaining S(U), from finite-size calculations, these re- 
sults should be viewed with some caution. Experimental 
searches for such a behavior would be quite interesting. 



netic phase at large U/t and a semi-metal phase at small 
U/t. Quantum monte carlo results for the compressibil- 
ity, which looks at the charge response of the system, and 
the magnetic structure factor, which looks at the spin re- 
sponse, both suggest a transition around U c /t w 5. The 
series expansion results for the sublattice magnetization, 
which is the spin order parameter and the charge excita- 
tion gap, which characterizes the Mott transition, both 
point to a single transition at U/t s=s 4. The discrepency 
between the quantum monte carlo and series expansion 
results reflects the uncertainties in the calculations, espe- 
cially as the critical point is approached. Thus we expect 
the transition to lie in the range 4 < U/t < 5, a result in 
complete agreement with the previous work of Martelo 
et al. @. 

Finally, one of the goals of this work was to look for 
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TABLE II: Series coefficients for the spin-wave excitation spectrum A(k x , k y )/U and 1-hole dispersion Aih(fca;, k y )/U. Nonzero 
coefficients up to order A 13 for t/U — 0.15 and J/U — 0.0225 are listed. 
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